Introdução à Análise Estatística Espacial

O que é Análise Estatística Espacial ?

  • São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;

  • Defini-se “Análise estatística espacial quando os dados são espacialmente localizados e se considera explicitamente a possível importância de seu arranjo espacial na análise ou interpretação dos resultados” (Bailey & Gatrell, 1995).

  • A primeira pergunta a ser feita é: A distribuição dos dados apresenta um padrão aleatório ou apresenta uma agregação definida (clusters) ?

Origem da Estatística Espacial

  • Dr. John Snow (1813-1858) - Considerado pai da Epidemiologia Moderna

Mapeamento dos casos de coléra (\(\bullet\)) e as bombas de água (X) em Londres, 1854.

Objetivos da Estatística Espacial

  1. Análise de padrões espaciais e espaço-temporais (ex: Análise Exploratória Espacial, Correlação Espacial)

  2. Modelagem Espacial (ex: Modelos Estatísticos de Regressão Espacial)

Dependência Espacial ou Autocorrelação Espacial

  • “Independência é um pressuposto muito conveniente que faz grande parte da teoria estatística matemática tratável. Entretanto, modelos que envolvem dependência estatística são frequentemente mais realísticos. . . . Neste tipos de dados a dependência está presente em todas as direções e fica mais fraca a medida em que aumenta a dispersão na localização dos dados.” (Cressie N.,1991)

  • “Todas as coisas se parecem, coisas mais próxmas são mais parecidas que aquelas mais distantes.” (Tobler, 1979). Também conhecida como \(1^a\) Lei da Geografia

Tipologia dos Dados Espaciais

Os diferentes tipos de dados espaciais são tradicionalmente classificados de acordo com uma tipologia. Esta caracterização diz respeito a natureza estocástica da observação.

Cressie divide a estatı́stica espacial em 3 grandes áreas:

  1. Dados de processos pontuais

  2. Dados de geoestatística

  3. Dados de área

Existem métodos estatı́sticos diferentes para descrever ou analisar estes tipos de dados.

Padrões Pontuais

  • O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.

  • Exemplos: Localização de crimes, localização da residência dos casos de dengue, localização de espécies vegetais, etc.

  • Neste caso, o dado aleatório de interesse é a localização espacial do evento.

Estimativas de Kernel (ou Mapas de Calor)

Localização da ocorrência de casos de dengue em Belo Horizonte/MG

Geoestatística

  • Podemos definir como sendo uma análise de um atributo espacialmente contínuo amostrado em localizações fixas.

  • Os dados compreendem um conjunto de localizações (em geral latitudes e longitudes), mas agregados a eles uma medida, como por exemplo:

  • Bastante utilizada em ciências ambientais (ex: volume de chuva, concentração de poluentes no ar, número de ovos de Aedes postado em ovitrampas, etc.)

Exemplo: Mapa sobre o teor de argila no solo.

  • Aqui os pesos são determinados a partir de uma análise espacial, baseada no semivariograma experimental.

  • Além disso, a Krigeagem fornece, em média, estimativas não tendênciosas e com variância mínima.

Dados de Área

  • Na análise de áreas o atributo estudado é em geral resultado de uma contagem ou uma medida de síntese (em geral somas ou médias).

  • O objetivo é a detecção e explicação de padrões e tendências observados entre áreas.

  • Área é definida por um polígono cuja forma pode ser complexa bem como as relações de vizinhança.

Mapa Temático

Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998

Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998

Matriz de Vizinhança

Lisa Map ou Moran Local

Mapa temáticos com o Moran Global da Incidência de Hansenı́ase São Paulo,1989

Mapa temáticos com o Moran Global da Incidência de Hansenı́ase São Paulo,1989

Lisa Map ou Moran Local

Modelos de Regressão Espacial

  • Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial

  • Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);

  • Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).

  • Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:

    1. Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR): atribuem a autocorrelação ao erro;

    2. Locais: parâmetros variam continuamente no espaço (ex: Geographically weighted regression (GWR))

Exemplo com os dados de dengue em Dourados/MS

Nesta apresentação serão utilizados os dados que fizeram parte da dissertação de Isis Rodrigues Reitman intitulada “SAÚDE E AMBIENTE URBANO: A RELAÇÃO DE INCIDÊNCIA DE DENGUE E AS DISPARIDADES ESPACIAIS EM DOURADOS – MS”, apresentada no Programa de Pós-Graduação em Geografia da Universidade Federal da Grande Douradosos/MS, em abril de 2016.

Essa aplicação também se encontra no Curso de Estudos Ecológicos ministrado para o curso de Pós-Graduação em Epidemiologia em Saúde Pública em 2019, pelos pesquisadores Oswaldo Gonçalves Cruz (PROCC/FIOCRUZ) e Wagner Tassinari (DEMAT/ICE/UFRRJ)

Biliotecas do R que serão utilizadas

library(readr)
library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)

Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS

unzip("dados/pop2010.zip", exdir = "dados")
pop2010 <- read_csv("dados/pop2010.csv")

unzip("malhas/Setor_UTM_SIRGAS.zip", exdir = "malhas")
setor.sf <- read_sf("malhas/Setor_UTM_SIRGAS.shp", crs = 31981)

unzip("malhas/contorno.zip", exdir = "malhas")
contorno.sf <- read_sf("malhas/contorno.shp", crs = 31981)

Fazendo um join com as tabelas com os setores censitários + popoulação

setor.sf <- setor.sf %>% mutate(idsetor = as.numeric(CD_GEOCODI)) %>% inner_join(pop2010, 
    by = "idsetor")

Lendo e plotando os casos de dengue georreferenciados em Dourados/MS

unzip("dados/dengue_dourados.zip", exdir = "dados")
casos <- read_csv("dados/dengue_dourados.csv")
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = casos.sf, 
    color = "red", size = 1) + theme_void()

Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS

unzip("dados/lixo_dourados.zip", exdir = "dados")
lixo <- read_csv2("dados/lixo_dourados.csv")
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo.sf, 
    color = "blue", size = 1) + theme_void()

Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS

Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.

lixo2.sf <- st_intersection(contorno.sf, lixo.sf)
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo2.sf, 
    color = "blue", size = 1) + theme_void()

Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf, color = "red", 
    size = 0.7, aes(colour = "Caso"), show.legend = "point") + geom_sf(data = lixo2.sf, 
    color = "salmon", size = 1, aes(colour = "Lixo"), show.legend = "point") + scale_fill_distiller(palette = "PuBu", 
    direction = 1) + scale_colour_manual(values = c(Caso = "red", Lixo = "slamon")) + 
    theme_minimal()

Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500)

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = lixo_buffer, color = "gray", 
    fill = "transparent", size = 0.4) + geom_sf(data = casos.sf, color = "red", size = 0.7) + 
    scale_fill_distiller(palette = "PuBu", direction = 1) + scale_colour_manual(values = c(Caso = "red", 
    Lixo = "slamon")) + theme_minimal()

Represntando os casos e o lixo de forma interativa

tm_shape(setor.sf) + tm_borders("black") + tm_shape(casos.sf) + tm_dots("red") + 
    tm_shape(lixo.sf) + tm_dots("green") + tm_shape(lixo_buffer) + tm_borders("blue") + 
    tmap_mode("view")

Convertendo o dado de pontos (padrão pontual) para dados de área

## conta casos por setor 
casos.sf$contador <- 1 
casos <- setor.sf %>%  
  st_join(casos.sf) %>% 
  filter(CLASSI_FIN == 1) %>%  ## seleciona somente os casos confirmados
  group_by(ID) %>% 
  summarise(casos = sum(contador))

st_geometry(casos) <- NULL  ## remove atributos de geometria

## numero de depositos de Lixo por setor 
lixo.sf$contador <- 1 

lixo <- setor.sf %>%  
  st_join(lixo.sf) %>% 
  group_by(ID) %>% 
  summarise(lixo = sum(contador))

st_geometry(lixo) <- NULL ## remove atributos de geometria

# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf %>% 
  left_join(lixo,by = 'ID') %>% 
  left_join(casos,by = 'ID') 

Plotando o mapa temático dos casos por setor censitário

plot(setor.sf["casos"])

Plotando o mapa temático dos pontos de coleta de lixo por setor censitário

plot(setor.sf["lixo"])

Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário

setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0  # Transformando os missings em zero

summary(setor.sf$tx)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.736   4.065   4.311  56.399 

Plotando a distribuição da incidência em Dourados/MS

library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + geom_sf(aes(fill = tx), color = "black") + scale_fill_gradientn(colours = pal) + 
    ggtitle("Taxa de incidência de Dengue") + theme_void()

Kernel por atributos

Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.

Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

dourados.w <- as.owin(st_geometry(dourados.contorno))

Extraindo os centróidas dos polígonos em Dourados/MS

centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")

plot(centroides.sp)

Colocando os pontos no formato sp

centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y, dourados.w)

plot(centroides.ppp, pch = 19, cex = 0.5)

Fazendo o kernel por atributo da taxa de detecção

kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)

Matriz de Vizinhança

Construindo a matriz de vizinhança para verificar a autocorrelação espacial

library(spdep)
viz <- poly2nb(setor.sf)
viz
Neighbour list object:
Number of regions: 284 
Number of nonzero links: 1726 
Percentage nonzero weights: 2.139952 
Average number of links: 6.077465 

Iremos precisar da coordenadas dos centróides

setor.sp <- as(setor.sf, "Spatial")  # convertendo em formato sp
coord <- coordinates(setor.sp)  # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

Verificando a malha de conectividade da vizinhança de Dourados/MS

viz.sf <- as(nb2lines(viz, coords = coord), "sf")
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))

# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + geom_sf(fill = "salmon", color = "white") + geom_sf(data = viz.sf) + 
    theme_minimal() + ggtitle("Vizinhança por \n conectividade") + ylab("Latitude") + 
    xlab("Longitude")
mapa.viz

Obtendo a correlação da taxa de incidência de dengue Dourados/MS

pesos.viz <- nb2listw(viz)
moran.test(setor.sf$tx, pesos.viz)

    Moran I test under randomisation

data:  setor.sf$tx  
weights: pesos.viz    

Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.524720303      -0.003533569       0.001128660 

Plotando o correlograma

correl <- sp.correlogram(viz, setor.sf$tx, order = 8, method = "I")
correl
Spatial correlogram for setor.sf$tx 
method: Moran's I
           estimate expectation    variance standard deviate Pr(I) two sided
1 (284)  0.52472030 -0.00353357  0.00112866          15.7239       < 2.2e-16
2 (284)  0.23776816 -0.00353357  0.00048540          10.9525       < 2.2e-16
3 (284)  0.09536593 -0.00353357  0.00028358           5.8729       4.282e-09
4 (284) -0.02063378 -0.00353357  0.00019736          -1.2172         0.22351
5 (284) -0.07589158 -0.00353357  0.00016732          -5.5939       2.221e-08
6 (284) -0.06448448 -0.00353357  0.00016165          -4.7940       1.635e-06
7 (284) -0.02708340 -0.00353357  0.00016725          -1.8210         0.06861
8 (284) -0.02744543 -0.00353357  0.00018328          -1.7662         0.07735
           
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)    
5 (284) ***
6 (284) ***
7 (284) .  
8 (284) .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(correl)

Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.

setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[, 5]

tm_shape(setor.sf) + tm_polygons(col = "pval", title = "p-valores", breaks = c(0, 
    0.01, 0.05, 0.1, 1), border.col = "white", palette = "-Oranges") + tm_scale_bar(width = 0.15)

Moran Local (Lisa Map) da taxa de incidência Dourados/MS

resI <- localmoran.sad(lm(setor.sf$tx ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10, ]
    Local Morans I Stand. dev. (N)   Pr. (N) Saddlepoint  Pr. (Sad)
1 1   -0.009885292     -0.01575255 0.5062841 -0.03600714 0.51436167
2 2    0.226981101      0.46511580 0.3209243  0.70056077 0.24178858
3 3   -0.010910944     -0.01829621 0.5072987 -0.04041673 0.51611955
4 4    0.484675519      1.31014141 0.0950740  1.43930439 0.07503215
5 5    0.018648578      0.05952726 0.4762661  0.09384037 0.46261798
     Expectation  Variance    Skewness Kurtosis   Minimum  Maximum
1 1 -0.003533569 0.1625854 -0.05147857 8.753481 -57.75455 56.75455
2 2 -0.003533569 0.2456263 -0.04188294 8.752890 -70.87400 69.87400
3 3 -0.003533569 0.1625854 -0.05147857 8.753481 -57.75455 56.75455
4 4 -0.003533569 0.1388594 -0.05570264 8.753780 -53.41199 52.41199
5 5 -0.003533569 0.1388594 -0.05570264 8.753780 -53.41199 52.41199
            omega       sad.r       sad.u
1 1 -0.0001369880 -0.01569409 -0.01569909
2 2  0.0028119325  0.44472643  0.49831660
3 3 -0.0001590844 -0.01822753 -0.01823490
4 4  0.0066304485  1.08297547  1.59298211
5 5  0.0005595065  0.05929966  0.05942125
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]
setor.sf$MoranLocal <- summary(resI)[, 1]

library(scales)

ggplot(setor.sf) + geom_sf(aes(fill = MoranLocal), color = "black") + scale_fill_gradientn(colours = c("blue", 
    "white", "red"), values = rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))), 
    guide = "colorbar") + ggtitle("Moran local") + theme_void()

Modelos Linear, CAR e GWR

Ajustando o modelo de regressão linear simples.

setor.sf$lixo[is.na(setor.sf$lixo)] <- 0  # Transformando os missings em zero

dourados.lm <- lm(tx ~ lixo, data = setor.sf)
summary(dourados.lm)

Call:
lm(formula = tx ~ lixo, data = setor.sf)

Residuals:
   Min     1Q Median     3Q    Max 
-4.510 -4.115 -2.286  0.237 51.889 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.5103     0.4889   9.225   <2e-16 ***
lixo         -0.3955     0.1945  -2.033    0.043 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.365 on 282 degrees of freedom
Multiple R-squared:  0.01445,   Adjusted R-squared:  0.01095 
F-statistic: 4.134 on 1 and 282 DF,  p-value: 0.04298

Checando os residuos para verificar a presença de autocorrelação.

dourados.lm$lmresid <- residuals(dourados.lm)
moran.test(dourados.lm$lmresid, pesos.viz)

    Moran I test under randomisation

data:  dourados.lm$lmresid  
weights: pesos.viz    

Moran I statistic standard deviate = 15.578, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.520064269      -0.003533569       0.001129669 

Ajustando o modelo CAR

dourados.car <- errorsarlm(tx ~ lixo, data = setor.sf, listw = pesos.viz)
summary(dourados.car)

Call:errorsarlm(formula = tx ~ lixo, data = setor.sf, listw = pesos.viz)

Residuals:
      Min        1Q    Median        3Q       Max 
-12.53213  -2.34032  -1.09686   0.82021  31.62992 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  4.07199    1.54998  2.6271 0.008611
lixo        -0.25129    0.14814 -1.6963 0.089826

Lambda: 0.8063, LR test value: 167.98, p-value: < 2.22e-16
Asymptotic standard error: 0.041575
    z-value: 19.394, p-value: < 2.22e-16
Wald statistic: 376.13, p-value: < 2.22e-16

Log likelihood: -885.0643 for error model
ML residual variance (sigma squared): 25.435, (sigma: 5.0433)
Number of observations: 284 
Number of parameters estimated: 4 
AIC: 1778.1, (AIC for lm: 1944.1)

Checando os residuos para verificar a presença de autocorrelação

dourados.car$carresid <- residuals(dourados.car)
moran.test(dourados.car$carresid, pesos.viz)

    Moran I test under randomisation

data:  dourados.car$carresid  
weights: pesos.viz    

Moran I statistic standard deviate = 0.78357, p-value = 0.2166
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.023005067      -0.003533569       0.001147099 

Ajustando o modelo GWR

# Biblioteca para ajustar o modelos GWR
library(spgwr)
# Estimando a largura de banda “ideal” para o kernel
GWRbanda <- gwr.sel(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y), 
    adapt = T)
Adaptive q: 0.381966 CV score: 14699.46 
Adaptive q: 0.618034 CV score: 15071.4 
Adaptive q: 0.236068 CV score: 14204.85 
Adaptive q: 0.145898 CV score: 13336.85 
Adaptive q: 0.09016994 CV score: 12151.44 
Adaptive q: 0.05572809 CV score: 10995.75 
Adaptive q: 0.03444185 CV score: 10093.37 
Adaptive q: 0.02128624 CV score: 9758.837 
Adaptive q: 0.01315562 CV score: 9295.706 
Adaptive q: 0.008130619 CV score: 8996.895 
Adaptive q: 0.005024999 CV score: 8988.675 
Adaptive q: 0.00638842 CV score: 9000.179 
Adaptive q: 0.00310562 CV score: 9842.835 
Adaptive q: 0.004291861 CV score: 9067.819 
Adaptive q: 0.005545779 CV score: 8976.568 
Adaptive q: 0.005867639 CV score: 8980.638 
Adaptive q: 0.005586469 CV score: 8976.67 
Adaptive q: 0.005505089 CV score: 8976.598 
Adaptive q: 0.005545779 CV score: 8976.568 
dourados.gwr = gwr(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y), 
    adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)

dourados.gwr
Call:
gwr(formula = tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, 
    centroides.sp$Y), adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
Kernel function: gwr.Gauss 
Adaptive quantile: 0.005545779 (about 1 of 284 data points)
Summary of GWR coefficient estimates at data points:
                   Min.    1st Qu.     Median    3rd Qu.       Max.  Global
X.Intercept.  -0.324952   1.412188   3.099163   5.397652  37.581278  4.5103
lixo         -18.321709  -1.222217  -0.400929  -0.061492   1.555744 -0.3955
Number of data points: 284 
Effective number of parameters (residual: 2traceS - traceS'S): 134.4484 
Effective degrees of freedom (residual: 2traceS - traceS'S): 149.5516 
Sigma (residual: 2traceS - traceS'S): 5.447405 
Effective number of parameters (model: traceS): 100.4864 
Effective degrees of freedom (model: traceS): 183.5136 
Sigma (model: traceS): 4.917575 
Sigma (ML): 3.952993 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 1904.233 
AIC (GWR p. 96, eq. 4.22): 1687.144 
Residual sum of squares: 4437.827 
Quasi-global R2: 0.7140881 
# Colocando a saída do modelo dentro de um dataframe.

results <- as.data.frame(dourados.gwr$SDF)
head(results)
     sum.w X.Intercept.       lixo X.Intercept._se  lixo_se       gwr.e
1 4.630090     6.489949 -1.8954895        2.411594 1.718037 -0.02159391
2 3.934954     7.883872 -1.2811203        2.504116 1.630176  1.98454873
3 3.062602     7.267798 -1.5625806        2.780695 1.625888 -1.85906371
4 4.139484     9.791258 -1.1687314        2.376803 1.277302  6.80957181
5 5.541437     4.586352 -0.6839813        2.090863 1.185501 -2.63010618
      pred  pred.se   localR2 X.Intercept._se_EDF lixo_se_EDF pred.se.1
1 2.698970 2.477284 0.6403656            2.671425    1.903141  2.744191
2 7.883872 2.504116 0.5043081            2.773914    1.805815  2.773914
3 5.705218 2.164469 0.5288378            3.080292    1.801064  2.397673
4 8.622527 1.794393 0.4101644            2.632885    1.414921  1.987725
5 3.902371 1.524868 0.6665163            2.316137    1.313229  1.689160
   coord.x coord.y
1 731406.5 7541547
2 730845.4 7541481
3 730957.7 7541239
4 730487.2 7541437
5 729874.1 7541446
 [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

Verificando a distribuição dos coeficientes de regressão para a variável lixo

hist(results$lixo)
abline(v = median(results$lixo), col = "red")

Verificando a distribuição dos localR2

hist(results$localR2)
abline(v = median(results$localR2), col = "blue")

Incorporando alguns parâmetros de saída do modelo na tabela setor.sf

setor.sf$coef.lixo <- results$lixo
setor.sf$localR2 <- results$localR2
setor.sf$pred.gwr <- results$pred

Mapa dos coeficientes de regressão para a variável lixo

map.lixo <- ggplot(setor.sf) + geom_sf(aes(fill = coef.lixo), color = "black") + 
    scale_fill_gradientn(colours = pal) + ggtitle("Distribuição dos coef. var. lixo") + 
    theme_void()
map.lixo

Checando os residuos para verificar a presença de autocorrelação para o modelo GWR.

# Calculando os resíduos para o modelo GWR
results$residuos <- setor.sf$tx - results$pred

moran.test(results$residuos, pesos.viz)

    Moran I test under randomisation

data:  results$residuos  
weights: pesos.viz    

Moran I statistic standard deviate = 1.4806, p-value = 0.06935
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.046816835      -0.003533569       0.001156432 

Mapeando os coeficientes de regressão para a variável lixo por significancia através do teste de wald

# Calculando a estatística de wald
setor.sf$wald.teste <- abs(results$lixo/results$lixo_se)
# Dicotomizando a estatística de wald
setor.sf$wald.teste <- ifelse(setor.sf$wald.teste < 2, 0, 1)


map.wald <- ggplot(setor.sf) + geom_sf(aes(fill = factor(wald.teste)), color = "black") + 
    scale_fill_manual(values = c("white", "purple"), labels = c("< 2", ">= 2"), name = "Wald") + 
    ggtitle("Coef. lixo significativos") + theme_void()


library(gridExtra)
grid.arrange(map.lixo, map.wald, ncol = 2)

Mapa dos coeficientes de determinação regionalizados (\(R^2\) local).

ggplot(setor.sf) + geom_sf(aes(fill = localR2), color = "black") + scale_fill_gradientn(colours = pal) + 
    ggtitle("R² local") + theme_void()

Verificando a distribuição dos preditos.

histdens <- function(x, titulo = "") {
    densi <- density(x)
    xli <- range(densi$x)
    yli <- range(densi$y)
    hist(x, col = "red", probability = T, xlim = xli, ylim = yli, main = titulo)
    lines(densi, lwd = 2)
    abline(v = median(x), lwd = 4, col = 4, lty = 2)
}

attach(mtcars)
par(mfrow = c(2, 2))

hist.tx <- histdens(setor.sf$tx, "Tx Bruta")
hist.lm <- histdens(dourados.lm$fitted.values, "Pred LM")
hist.car <- histdens(dourados.car$fitted.values, "Pred CAR")
hist.gwr <- histdens(results$pred, "Pred GWR")

Mapeando os valores observados e preditos dos modelos ajustados

library(colorspace)  # 

setor.sf$brks <- cut(setor.sf$tx, include.lowest = TRUE, right = TRUE, breaks = c(-4, 
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))

tx.bruta <- ggplot(setor.sf) + geom_sf(aes(fill = brks), color = "black") + ggtitle("Taxa Bruta") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()



setor.sf$brks.lm <- cut(dourados.lm$fitted.values, lowest = TRUE, right = TRUE, breaks = c(-4, 
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.lm), color = "black") + ggtitle("Taxa Predita - LM") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()

setor.sf$brks.car <- cut(dourados.car$fitted.values, lowest = TRUE, right = TRUE, 
    breaks = c(-4, 0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", 
        "> 10"))


pred.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - CAR") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()

setor.sf$brks.gwr <- cut(results$pred, lowest = TRUE, right = TRUE, breaks = c(-4, 
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - GWR") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()

library(gridExtra)
grid.arrange(tx.bruta, pred.lm, pred.car, pred.gwr, ncol = 2)

Verificando a distribuição dos resíduoas .

library(vioplot)
vioplot(dourados.lm$residuals, dourados.car$residuals, results$residuos, names = c("LM", 
    "CAR", "GWR"), col = "orange")
title("Gráficos violinos da distribuição dos resíduos")

Mapeando os resíduos dos modelos ajustados

library(colorspace)  # 

setor.sf$brks.res.lm <- cut(dourados.lm$residuals, include.lowest = TRUE, right = TRUE, 
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5", 
        "> 5"))

res.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.lm), color = "black") + 
    ggtitle("Resíduos - LM") + scale_fill_discrete_sequential(palette = "Purple-Yellow", 
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", 
    drop = FALSE, name = "Taxa") + theme_void()

setor.sf$brks.res.car <- cut(dourados.car$residuals, include.lowest = TRUE, right = TRUE, 
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5", 
        "> 5"))

res.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.car), color = "black") + 
    ggtitle("Resíduos - CAR") + scale_fill_discrete_sequential(palette = "Purple-Yellow", 
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", 
    drop = FALSE, name = "Taxa") + theme_void()


setor.sf$brks.res.gwr <- cut(results$residuos, include.lowest = TRUE, right = TRUE, 
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5", 
        "> 5"))

res.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.gwr), color = "black") + 
    ggtitle("Resíduos - GWR") + scale_fill_discrete_sequential(palette = "Purple-Yellow", 
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", 
    drop = FALSE, name = "Taxa") + theme_void()


library(gridExtra)
grid.arrange(res.lm, res.car, res.gwr, ncol = 2)

Bibliografia sugerida

  • Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) . Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.

  • Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995

  • Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004

  • Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013

  • Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.